/*
 * Multiphase.cpp
 *
 *  Created on: 23 Jul 2011
 *      Author: Allan
 */

#include "Multiphase.h"

// C++ includes
#include <algorithm>
#include <iomanip>

// Boost includes
#include <boost/foreach.hpp>

Multiphase::Multiphase()
{}

Multiphase::Multiphase(const vector<Phase>& phases) :
phases(phases)
{
	// Loop over all the phases to set the sequence of species 
	BOOST_FOREACH(const Phase& phase, phases)
		species.insert(species.end(), phase.GetSpecies().begin(), phase.GetSpecies().end());
	
	species_map = SpeciesMap(species);
	
	// Initialiaze the species-phase mapping
	for(unsigned phase_index = 0; phase_index < phases.size(); ++phase_index)
		BOOST_FOREACH(const string& s, phases[phase_index].GetSpecies())
			species_to_phase[s] = phase_index;
}

const unsigned
Multiphase::GetNumPhases() const
{
	return phases.size();
}

const unsigned
Multiphase::GetNumSpecies() const
{
	return species.size();
}
	
const vector<Phase>&
Multiphase::GetPhases() const
{
	return phases;
}

const vector<string>&
Multiphase::GetSpecies() const
{
	return species;
}

const Index
Multiphase::FindPhase(const string& phase_name) const
{
	auto check_name_func = [&](const Phase& phase) { return phase.GetName() == phase_name; };
	
	return find_if(phases.begin(), phases.end(), check_name_func) - phases.begin();
}

const Index
Multiphase::operator[](const string& species) const
{
	return species_map[species];
}

const vector<Index>
Multiphase::operator[](const vector<string>& species) const
{
	return species_map[species];
}

const Index
Multiphase::PhaseIndex(const string& species) const
{
	return species_to_phase.find(species)->second;
}

template<class PhaseFunction> const VectorXd
GenericPhaseComputation(double T, double P, const VectorXd& n, const vector<Phase>& phases, const PhaseFunction& phase_function)
{
	// The activity-like data of the species to be computed
	VectorXd a(n.rows());
	
	// The current number of processed species, to be update at the end of every loop (see below)
	unsigned num_processed = 0;
	
	// Loop over all the phases
	BOOST_FOREACH(const Phase& phase, phases)
	{
		// The number of species in the current phase
		const unsigned num_species = phase.GetSpecies().size();
		
		// The vector of number of moles of the species of the current phase
		const VectorXd n_phase = n.segment(num_processed, num_species);
		
		// Compute the activity-like data of the species in the current phase using the template parameter phaseFunction
		const VectorXd a_phase = phase_function(phase, T, P, n_phase);
		
		// Transfer the activity-like data of the current phase to the vector of activity-like data a (beginning at position given by nprocessed)
		a.segment(num_processed, num_species) = a_phase;
		
		// Update the number of processed species
		num_processed += num_species;
	}

	return a;
}

const VectorXd
Multiphase::PhaseMoles(const Index& phase_index, const VectorXd& n) const
{
	unsigned offset = 0;
	
	for(unsigned pi = 0; pi < phase_index; ++pi)
		offset += phases[pi].GetNumSpecies();
	
	return n.segment(offset, phases[phase_index].GetNumSpecies());
}

const VectorXd
Multiphase::TotalMolesPhases(const VectorXd& n) const
{
	const unsigned num_phases = phases.size();
	
	VectorXd n_pi(num_phases);
	
	unsigned offset = 0;
	
	for(unsigned pi = 0; pi < num_phases; ++pi)
	{
		const unsigned num_species_phase_pi = phases[pi].GetNumSpecies();
		
		n_pi[pi] = n.segment(offset, num_species_phase_pi).sum();
		
		offset += num_species_phase_pi;
	}
	
	return n_pi;
}

const VectorXd
Multiphase::ActivityCoefficients(double T, double P, const VectorXd& n) const
{
	auto activity_coeffs_func = 
			[](const Phase& phase, double T, double P, const VectorXd& n) { return phase.ActivityCoefficients(T, P, n); };
	
	return GenericPhaseComputation(T, P, n, phases, activity_coeffs_func);
}
	
const VectorXd
Multiphase::Activities(double T, double P, const VectorXd& n) const
{
	auto activities_func = 
			[](const Phase& phase, double T, double P, const VectorXd& n) { return phase.Activities(T, P, n); };
	
	return GenericPhaseComputation(T, P, n, phases, activities_func);
}

const MatrixXd
Multiphase::ActivityDerivatives(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	const unsigned num_species = species.size();
	
	const unsigned num_phases = phases.size();
	
	MatrixXd derivatives = MatrixXd::Zero(num_species, num_species);
	
	unsigned offset = 0;
	
	for(unsigned iphase = 0; iphase < num_phases; ++iphase)
	{
		const unsigned num_species_phase = phases[iphase].GetNumSpecies();
		
		const VectorXd n_phase = n.segment(offset, num_species_phase);
		
		const VectorXd a_phase = a.segment(offset, num_species_phase);
		
		derivatives.block(offset, offset, num_species_phase, num_species_phase) = phases[iphase].ActivityDerivatives(T, P, n_phase, a_phase);
		
		offset += num_species_phase;
	}
	
	return derivatives;
}

ostream&
operator<<(ostream& out, const Multiphase& multiphase)
{
	BOOST_FOREACH(const Phase& phase, multiphase.GetPhases())
	{
		out << phase.GetName() << " Phase" << endl;
		
		BOOST_FOREACH(const string& s, phase.GetSpecies())
			out << "(" << multiphase[s] << ") " << s << endl;
		
		out << endl;
	}
	
	return out;
}
